/*
 * AqueousActivity.cpp
 *
 *  Created on: 19 Aug 2011
 *      Author: Allan
 */

#include "AqueousActivity.h"

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "AqueousModels/HkfModel.h"
#include "AqueousModels/SetschenowModel.h"
#include "AqueousModels/DrummondCO2Model.h"
#include "../ElectrolyteSolution.h"
#include "../Numerics/LinearInterpolator.h"
#include "../Utilities/Algorithms.h"

AqueousActivity::AqueousActivity()
{}

AqueousActivity::AqueousActivity(const vector<string>& species) :
AuxiliaryActivity(species)
{
	// Set the index of the water species
	iH2O = IndexOf("H2O", species);
	
	// The electrolyte solution used to create the activity coefficient functions
	ElectrolyteSolution sol(species);
	
	// The ions in the given list of (species)
	const vector<string> ions = sol.GetIons();
	
	// Initialize the activity coefficient functions of all the species with Setschenow's model (b = 0.1) and then change it for some species
	activitycoeffs.resize(species.size(), SetschenowModel(sol, 0.1));
	
	// Set the water Hkf model for H2O
	SetActivityCoefficient("H2O", HkfWaterModel(sol));
	
	// Set the ionic Hkf model for the ionic species
	BOOST_FOREACH(const string& ion, ions)
		SetActivityCoefficient(ion, HkfIonicModel(sol, ion));
	
	// Set the Drummond model for CO2(aq) (only effective if it was provided)
	SetActivityCoefficient("CO2(aq)", DrummondCO2Model(sol));
}

AqueousActivity::~AqueousActivity()
{}

const VectorXd AqueousActivity::Activities(double T, double P, const VectorXd& n) const
{
	VectorXd a(n.rows());
	
	// The molality factor to trasnform number of moles into molality
	const double molalityfactor = 55.508 / n[iH2O];
	
	for(int i = 0; i < a.rows(); ++i) 
		a[i] = molalityfactor * n[i] * activitycoeffs[i](T, P, n);
	
	return a;
}

const ActivityCoefficient IdealAqueousModel(const ElectrolyteSolution& sol)
{
	/// The index of the water species
	const Index iH2O = sol.GetWaterIndex();
	
	auto molarfractionWater = [=](double T, double P, const VectorXd& n) 
			{ return n[iH2O] / n.sum(); };
	
	return molarfractionWater;
}

const ActivityCoefficient HkfIonicModel(const ElectrolyteSolution& sol, const string& ion)
{
	// The ions in the given electrolyte solution (sol)
	const vector<string> ions = sol.GetIons();
	
	// The effective electrostatic radii of the ions
	vector<double> effRadii;
	
	BOOST_FOREACH(const string& ion, ions) 
		effRadii.push_back(EffectiveElectrostaticRadius(ion));
	
	// The index of the given (ion) among the (ions) in the electrolyte solution 
	const Index idxion = IndexOf(ion, ions);
	
	// Generate the activity coefficient function for the given ion
	return std::bind(HkfIonicActivityCoefficient, sol, _1, _2, _3, idxion, effRadii);
}

const ActivityCoefficient HkfWaterModel(const ElectrolyteSolution& sol)
{
	// The ions in the given electrolyte solution (sol)
	const vector<string> ions = sol.GetIons();
	
	// The effective electrostatic radii of the ions
	vector<double> effRadii;
	
	BOOST_FOREACH(const string& ion, ions)
		effRadii.push_back(EffectiveElectrostaticRadius(ion));
	
	// Generate the activity coefficient function for the water species
	return std::bind(HkfWaterActivityCoefficient, sol, _1, _2, _3, effRadii);
}

const ActivityCoefficient SetschenowModel(const ElectrolyteSolution& sol)
{
	return SetschenowModel(sol, 0.0);
}

const ActivityCoefficient SetschenowModel(const ElectrolyteSolution& sol, double b)
{
	return SetschenowModel(sol, {b}, {25.0}); // Assume the value b is given at 25.0 degC
}

const ActivityCoefficient SetschenowModel(const ElectrolyteSolution& sol, const vector<double>& bi, const vector<double>& Ti)
{
	LinearInterpolator b(Ti, bi);
	
	return std::bind(SetschenowActivityCoefficient, sol, _1, _2, _3, b);
}

const ActivityCoefficient DrummondCO2Model(const ElectrolyteSolution& sol)
{
	return std::bind(DrummondCO2ActivityCoefficient, sol, _1, _2, _3);
}

